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Using quantum Monte Carlo (QMC) simulations we study the ground-state properties of the one- 
dimensional fermionic Hubbard model in traps with an underlying lattice. Since due to the confining 
potential the density is space dependent, Mott-insulating domains always coexist with metallic 
regions, such that global quantities are not appropriate to describe the system. We define a local 
compressibility that characterizes the Mott-insulating regions and analyze other local quantities. 
It is shown that the momentum distribution function, a quantity that is commonly considered in 
experiments, fails in giving a clear signal of the Mott-insulator transition. Furthermore, we analyze 
a mean-field approach to these systems and compare it with the numerically exact QMC results. 
Finally, we determine a generic form for the phase diagram that allows us to predict the phases to 
be observed in the experiments. 
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I. INTRODUCTION 

The realization of Bose-Einstein condensation (BEC) 
of trapped atomic gases Q, has generated in the 

last years a huge amount of experimental and theoretical 
research 0, 1^. BEC was achieved by confining atoms 
in magnetic traps and lowering the temperature of the 
system via the evaporative cooling technique. In this 
technique, the hottest atoms are selectively removed from 
the system and the remaining ones rethermalize via two- 
body collisions. A common feature of these experiments 
is that the trapped gases are dilute; mean-field theory 
then provides a useful framework to study the role of the 
interaction between particles. 

Recently, a new feature has been added to the exper- 
iments: the magnetically trapped condensate is trans- 
ferred into an optical lattice generated by interfering 
laser beams 0. With this new experimental setup it 
is possible to access the strongly correlated regime. The 
superfluid-Mott-insulator transition was studied and 
other interesting physical phenomena such as the col- 
lapse and revival of the condensate have been found 0. 
The presence of the optical lattice and the fact that the 
particles interact only via contact interaction lead in a 
natural way to the Hubbard model as a paradigm for 
these systems. A theoretical work proposing such exper- 
iments M and recent quantum Monte Carlo (QMC) simu- 
lations [aEl have investigated these systems beyond the 
mean-field approximation. It has been found that the in- 
compressible Mott-insulating phase always coexists with 
compressible phases so that a local order parameter 
has to be defined to characterize the system. Also, the 
phase diagram of trapped bosons was found to be more 
complex than the one in the homogeneous system |^. 

The experimental realization of ultracold fermionic 
gases is more difficult than the bosonic case, and has 
been achieved only recently [H IH 111 El EE El 113 . 
Unlike bosons, single species fermions cannot be directly 
evaporatively cooled to very low temperatures because 
the s-wave scattering that could allow the gas to rether- 



malize during the evaporation is prohibited for identi- 
cal fermions. This problem has been overcome by si- 
multaneously trapping and evaporatively cooling two- 
component Fermi gases IT^ , and introducing mixed 
gases of bosons and fermions in which bosons enable 
fermions to rethermalize through their elastic interac- 
tions [l^ ^13,, . More recently, rapid forced evap- 
oration employing a Feshbach resonance of two different 
spin states of the same fermionic atom has been used 
[l7l |. Now that it is possible to go well below the degen- 
eracy temperature in the experiments and superfluidity 
appears within reach E3i it is expected that the metal- 
Mott-insulator transition (MMIT) could also be realized 
for fermions on an optical lattice. 

Motivated by this expectation we study the ground 
state of the one-dimensional (ID) fermionic Hubbard 
model with a harmonic trap and with repulsive contact 
interaction, using QMC simulations. Like in the bosonic 
case 111, Mott-insulating domains appear over a continu- 
ous range of fillings and always coexist with compressible 
phases, such that global quantities are not appropriate to 
characterize the system. Instead, we define a local-order 
parameter (a local compressibility) in order to charac- 
terize the local phases present in the system. We also 
analyze the generic features that are valid for any kind 
of confining potential, and not only for the harmonic one. 

It is well known that the Hubbard model in the peri- 
odic case displays a MMIT phase transition at half filling 
and at a finite value of the on-site repulsive interaction 
U. (In the case of perfect nesting the transition occurs 
at t/ = 0.) The two routes to this MMIT are the fiUing- 
controlled MMIT and the bandwidth-controlled MMIT 
El|. In the confined case, it is possible to drive the tran- 
sitions by changing the total filling of the trap, the on- 
site repulsive interaction, or varying the curvature of the 
confining potential. In the following section, we study 
a number of local quantities as a function of the fill- 
ing and the strength of the interaction. Although they 
signal the appearance of a Mott-insulating phase, such 
quantities provide no rigorous criterion to characterize 
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the Mott-insulating region. We therefore define a local 
compressibility that acts as a genuine local-order param- 
eter. In Sec. Ill, we study the momentum distribution 
function for the confined system. In Sec. IV, we discuss 
a mean-field approach for the ID trapped system and 
compare the results with the ones obtained using QMC. 
In Sect. V, we analyze the phase diagram and determine 
its generic form, which allows to compare systems with 
different sizes, number of particles and curvatures of the 
harmonic confining potential. In this section we also an- 
alyze the extension, for arbitrary confining potentials, of 
the results obtained for the harmonic case. Finally, the 
conclusions are given in Sec. VI. Some of the results pre- 
sented here were summarized in a previous publication 

II. LOCAL ORDER PARAMETER 

Due to the fact that the confining potential leads to an 
inhomogeneous density profile, we study in the present 
section how local quantities behave in the trapped system 
when the parameters at hand are changed. The Hamil- 
tonian of the fermionic Hubbard model with a confining 
parabolic potential has the form 

i.(T i 

+Vj2^'n,,, (1) 

where cj^, c^^ are creation and annihilation operators, 
respectively, for a fermion with spin a at site i, and 
riia- — cl^Cia^ such that t is the hopping amplitude, U 
is the on-site interaction that in the present work will 
be considered repulsive {U > 0), V is the curvature 
of the confining harmonic potential, and Xi measures 
the position of the site i {xi = ia with a the lattice 
constant). The number of lattice sites is N and is se- 
lected so that all the fermions are confined in the trap. 
We denote the total number of fermions in the trap as 
Nf and consider equal number of fermions with spins 
up and down {Nf^ = Nfi = Nf/2). In our simula- 
tions, we used the zero-temperature projector method 
|20l l2l| adapted from the QMC determinantal algorithm 
by Blankenbecler, Scalapino, and Sugar [2^l2^l2J |. The 
discrete Hubbard-Stratonovich transformation by Hirsch 
was used [2^ . For details of the algorithm we refer to a 
number of reviews [2fiLl27j. 

The results for the evolution of the local density {rii ~ 
{nil + ^ function of the total number of the 

confined particles, are shown in Fig. ^ For the lowest 
filling, so that n < 1 at every site, the density shows a 
profile with the shape of an inverted parabola, similar to 
that obtained in the non-interacting case |28l |. and hence, 
such a situation should correspond to a metallic phase. 
Increasing the number of fermions up to Nj = 60, a 
plateau with n ~ I appears in the middle of the trap. 



surrounded by a region with n < 1 (metallic). Since in 
the homogeneous case, a Mott insulator appears at n = 1, 
it is natural to identify the plateau with such a phase. 
The Mott-insulating domain in the center of the trap 
increases its size when more particles are added, but at a 
certain filling {Nf ~ 70 here) this becomes energetically 
unfavorable and a new metallic phase with n > 1 starts to 
develop in the center of the system. Upon adding more 
fermions, this new metallic phase widens spatially and 
the Mott-insulating domains surrounding it are pushed to 
the borders. Depending on the on-site repulsion strength, 
they can disappear and a complete metallic phase can 
appear in the system. Finally, a "band insulator" (i.e., 
n = 2) forms in the middle of the trap for the highest 
fillings (after Nf = 144 here). Due to the full occupancy 
of the sites, it will widen spatially and push the other 
phases present in the system to the edges of the trap 
when more fermions are added. 




FIG. 1: Evolution of the local density in a parabolic con- 
fining potential as a function of the position in the trap and 
increasing total number of fermions. The parameters involved 
are = 150, U ^ 4t and Va'^ = 0.002f. The positions are 
measured in units of the lattice constant a. 

Although the existence of flat regions in the den- 
sity profile is an indication that there is an insulator 
there, a more quantitative characterization is needed. 
As shown in Ref. |^, the variance of the local density 
{Ai ~ {nf) — {ni)^) may be used on a first approach 
(from here on, we refer to the variance as the variance of 
the local density) . In Fig. |2 we show four characteristic 
profiles present in Fig. ^ and their respective variances 
when the number of fermions in the system are Nf =50 
(a), 68 (b), 94 (c), and 150 (d). For the case in which 
only a metallic phase is present [Fig. |5fa)], it is possible 
to see that the variance decreases when the density ap- 
proaches n = 1 and has a minimum for densities close 
to that value. For the Mott-insulating domain [n = 1 
in Fig. Efb)], the variance has a constant value that is 
smaller than that of the metal surrounding it. As soon 
as the Mott-insulating phase is destroyed in the middle 
of the trap and a new metallic region with n > 1 devel- 
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ops there [Fig.[2Ic)], the variance increases in this region. 
The variance in the metaUic region with n > 1 will start 
to decrease again when the density approaches n = 2 
and will have values even smaller than those in the Mott- 
insulating phase. Finally, when the insulator with n ~ 2 
is formed in the center of the trap the variance vanishes 
there [Fig.l^d)]. 
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FIG. 2: (Color online) Four density profiles (A) (cuts across 
Fig-0 £ind their variances (O)- The fillings are Nf —50 (a), 
68 (b), 94 (c), and 150 (d). 

Alternatively, Mott-insulating regions can be obtained 
by increasing the ratio between the on-site repulsive in- 
teraction and the hopping parameter as was done in 
experiments for bosons confined in optical lattices 0. 
Fig. Ola) shows the evolution of the density profiles 
in a trapped system with N = 100, Nf = 70, and 
Va^ = 0.0025t when this ratio is increased from U/t = 2 
to U/t = 8 [for more details of these density profiles see 
Fig. [ZIc)-[3h)] . It can be seen that for small values of 
U/t {U/t — 2) there is only a metallic phase present in 
the trap. As the value of U/t {U/t — 4) is increased, a 
Mott-insulating phase tries to develop at n = 1 while a 
metallic phase with n > 1 is present in the center of the 
system. As the on-site repulsion is increased even further 
{U/t = 6, 8), a Mott-insulating domain appears in the 
middle of the trap suppressing the metallic phase that 
was present there. In Fig. IHIb) we show the variance of 
the density for the profiles in Fig. Ola) (from top to bot- 
tom, the values presented are for C//i ~ 2, 4, 6, 8). As 
expected, the variance decreases in both the metallic and 
Mott-insulating phases when the on-site repulsion is in- 
creased. When the Mott-insulating plateau is formed in 
the density profile, a plateau with constant variance ap- 
pears in the variance profile with a value that will vanish 
only in the limit U/t — > oo. As shown in Fig.|3Jb). when- 
ever a Mott-insulating domain is formed in the trap, the 
value of the variance in it is exactly the same as the one 



for the Mott-insulating phase in the homogeneous system 
for the same value of U/t (horizontal dashed lines). This 
would support the validity of the commonly used local 
density (Thomas- Fermi) approximation |29j. However, 
the insets in Fig. |3fb), show that this is not necessarily 
the case, since for U/t = 4, the value of the variance in 
the Mott-insulating phase of the homogeneous system is 
still not reached in the trap, although the density reaches 
the value r? = 1. Therefore, in contrast to the homoge- 
neous case, a Mott-insulating region is not determined by 
the filling only. In the cases oiU/t = 6 [inset in Fig. ^h) 
for a closer look] and U/t = 8, the value of the variance 
in the homogeneous system is reached and then we can 
say that Mott-insulating phases are formed there. 




FIG. 3: (Color online) Profiles for a trap with Va^ = 0.0025i 
and Nf = 70, the on-site repulsions are U/t =2 (A), 4 (v)i 
6 (O), and 8 (O)- (s-) Local density, (b) variance of the 
local density, (c) local compressibility ft* as defined in Eq. 

The dashed lines in (b) are the values of the variance in 
the n = 1 homogeneous system for U/t = 2, 4, 6, 8 (from 
top to bottom). 



Although the variance gives a first indication for the 
formation of a local Mott insulator, an ambiguity is 
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still present, since there are metallic regions with den- 
sities very close to n = and n = 2, where the vari- 
ance can have even smaller values than in the Mott- 
insulating phases. Therefore, an unambiguous quantity is 
still needed to characterize the Mott-insulating regions. 
We propose a local compressibility local-order pa- 
rameter to characterize the Mott-insulating regions, that 
is defined as 

^4= Xt,t+j , (2) 

where 

Xt,j = {niUj) - (rij) {uj) (3) 

is the density-density correlation function and t{U) ~ 
b^{U), with £,{U) the correlation length of Xi.j in the un- 
confined system at half-filling for the given value of U. 
As a consequence of the charge gap opened in the Mott- 
insulating phase at half filling in the homogeneous sys- 
tem, the density-density correlations decay exponentially 
[X{x) oc exp~'^/^('^)] enabhng ^([/) to be determined. The 
factor b is chosen within a range where becomes qual- 
itatively insensitive to its precise value. Since there is 
some degree of arbitrariness in the selection of 6, we will 
examine its role in the calculation of here. We show 
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FIG. 4: (Color online) Local compressibility as a function of ^ 
for the homoeeneous and trapped {Va'^ = 0. 00254, Nf = 70) 
systems in the Mott-insulating and metallic phases, U = 6t. 
(v) n — 1 homogeneous, (O) ni=o = 1 trapped, (A) n = 0.66 
homogeneous, and (O) ni=_32 ^ 0.66 trapped. 

in Fig. ^ (in a semi-log plot) how the local compress- 
ibility behaves as a function of i for a homogeneous and 
a trapped systems in the Mott-insulating (n = 1) and 
metallic (n ^ 0.66) phases. In the Mott-insulating phase, 
it is possible to see that in the homogeneous and con- 
fined systems the local compressibility decays exponen- 
tially to zero in exactly the same way due to the charge 
gap present there. In the metallic case there is no charge 
gap so that density-density correlations decay as a power 
law and the local compressibility, shown in Fig. 01 de- 
cays very slowly. In Fig. 01 a departure in the behavior 
of the local compressibility for the trapped system from 
the homogeneous case for large values of i can be seen. 
This occurs because correlations with points very close 



to the band insulating (n = 0) and the Mott-insulating 
(n = 1) regions surrounding the metallic phase [see Fig. 
|3Ja)] are included in k. However, this effect does not af- 
fect the fact that the value of the local compressibility is 
clearly different from zero as long as the size of the metal- 
lic phase is larger than 2£. We obtain that considering 
6 ~ 10 (no matter what the exact value of b is), the lo- 
cal compressibility is zero with exponential accuracy (i.e. 
^ 10~® or less) for the Mott-insulating phase and has a 
finite value in the metallic phase. For the case in Fig. 0| 
where U = 6t, the correlation length is ^ ^ 0.8 so that 
^ 8. Physically, the local compressibility defined here 
gives a measure of the change in the local density due to 
a constant shift of the potential over a finite range but 
over distances larger than the correlation length in the 
unconfincd system. 

In Fig.|3{c), we show the profiles of the local compress- 
ibility for the same parameters as Figs.l^^a) and (b) (we 
did not include the profile of the local compressibility for 
U = 2t because for that value of U we obtain that £ is 
bigger than the system size). In Fig. |3{c), it can be seen 
that the local compressibility only vanishes in the Mott- 
insulating domains. For U = At, it can be seen that in 
the region with n ~ 1 the local compressibility, although 
small, does not vanish. This is compatible with the fact 
that the variance is not equal to the value in the homoge- 
neous system there, so that although there is a shoulder 
in the density profile, this region is not a Mott insulator. 
Therefore, the local compressibility defined here serves as 
a genuine local order parameter to characterize the insu- 
lating regions that always coexist with metallic phases. 
At t his p oint we would like to remark that, as shown in 
Ref. [la, the local compressibility shows critical behav- 
ior on entering the Mott-insulating region. This is the 
reason to speak about phases, since on passing from one 
region to the other, critical behavior sets in. 

Finally, wc discuss in this section the spin-spin cor- 
relation function, since in periodic chains, quasi- long- 
range antiferromagnetic correlations appear in the Mott- 
insulating phase. In Fig. Owe show the local {S^Sj) cor- 
relation function for some points of the profiles presented 
in Fig. 12] We measured the spin-spin correlations in the 
trap at the points in the figures where it can be clearly 
seen that {S^Sj) has the maximum value. In Fig. [SJa) 
it can be seen that in the metallic phase with n < 1, the 
spin-spin correlations decay rapidly and do not show any 
clear modulation, which is due to the fact that the density 
is changing in this region. In the local Mott-insulating 
phase [center of the trap in Fig.|3^b)], short-range antifer- 
romagnetic correlations are present, and they disappear 
completely only on entering in the metallic regions. For 
the shoulders with n ~ 1 [Fig. [3c)], the antiferromag- 
netic correlations are still present but due to the small 
size of these regions they decay very rapidly. Finally for 
the metallic regions with n > 1, the spin correlations be- 
have like in the metallic phases with n < 1 as shown in 
Fig.Hd). 
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FIG. 5: (Color online) Local spin-spin correlations in a trap 
measured with respect to the points in which in the figure they 
have their maximum value. The density at each point (total 
filling in the trap) is 71^=25 = 0.61 (A^/ =50) (a), ni=o = 1-00 
(Nj =68) (b), ni=25 = 1.01 (iV/ =94) (c), and ni=25 = 1-44 
(iV/ =150) (d), for a trap with N = 150, U = 4i, and Va^ = 
0.002t (the density and variance profiles for these parameters 
were presented in Fig. I^J. 



III. MOMENTUM DISTRIBUTION FUNCTION 

In most experiments with quantum gases carried out 
so far, the momentum distribution function, which is de- 
termined in time-of-flight measurements, played a cen- 
tral role. A prominent example is given by the study of 
the superfluid-Mott-insulator transition in the bosonic 
case. Also a QMC study relating this quantity to the den- 
sity profiles and proposing how to determine the point 
at which the superfluid-Mott-insulator transition occurs 
was presented in Ref. ^3 ■ As shown below, we find that 
this quantity is not appropriate to characterize the phases 
of the system in the fermionic case, and does not show 
any clear signature of the MMIT. 

In Fig.|Slwe show the normalized momentum distribu- 
tion function (n^) for the same density profiles presented 
in Fig. 121 For the trapped systems, we always normal- 
ize the momentum distribution to be unity at fc = 0. In 
the case presented in Fig. |S1 transitions between phases 
occur due to changes in the total filling of the trap. We 
first notice that for the pure metallic phase in the har- 
monic trap [Fig.EIa)] does not display any sharp feature 
corresponding to a Fermi surface, in clear contrast to the 
homogeneous case. The lack of a sharp feature for the 
Fermi surface is independent of the presence of the inter- 
action and is also independent of the size of the system. 
In the non-interacting case, this can be easily understood: 
the spatial density and the momentum distribution will 



have the same functional form because the Hamiltonian 
is quadratic in both coordinate and momentum. When 
the interaction is present, it could be expected that the 
formation of local Mott-insulating domains generates a 
qualitatively and quantitatively different behavior of the 
momentum distribution, like in the homogeneous case 
where in the Mott-insulating phase the Fermi surface dis- 
appears and rifc is smoother. In Fig.|S{b), it can be seen 
that there is no qualitative change of the momentum dis- 
tribution when the Mott-insulating phase is present in 
the middle of the trap. Quantitatively in this case 
is very similar to the pure metallic cases Fig. Ela) (for 
densities smaller than one) and Fig. (HJc) (where in the 
middle of the trap the density is higher than one). Only 
when the insulator with n = 2 appears in the middle of 
the system we find a quantitative change of for U = At, 
as shown in Fig. El^d). 
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FIG. 6: (Color online) Normalized momentum distribution 
function for the same parameters of Fig. |21 The fillings are 
Nf =50 (a), 68 (b), 94 (c), and 150 (d). 



We also studied the momentum distribution function 
when the MMIT is driven by the change of the on-site 
repulsion. We did not observe any clear signature of the 
formation of the Mott-insulating phase in rifc. In Fig. 
13 we show the normalized momentum distribution func- 
tion [Fig.|3;a)-C|;d)] for density profiles [Fig.[Z|;c)-Cl;h)] in 
which the on-site repulsion was increased from U/t — 2 
to 8. It can be seen that the same behavior present in 
Fig. El and the quantitative changes in appear only 
when the on-site repulsion goes to the strong-coupling 
regime, but this is long after the Mott-insulating phase 
has appeared in the system. 

At this point one might think that in order to study 
the MMIT using the momentum distribution function, 
it is necessary to avoid the inhomogcneous trapping po- 
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FIG. 7: (Color online) Normalized momentum distribution 
function (a)-(d) and their corresponding density profiles (e)- 
(h) for U/t = 2 (a),(e), 4 (b),(f), 6 (c),(g), 8 (d),(h) and 
iV = 100, Nf = 70, Va^ = 0.0025t. 



tcntial and use instead a kind of magnetic box with in- 
finitely high potential on the boundaries. However, in 
that case one of the most important achievements of the 
inhomogeneous system is lost, i.e., the possibility of cre- 
ating Mott-insulating phases for a continuous range of 
fillings. In the perfect magnetic box, the Mott-insulating 
phase would only be possible at half filling, which would 
be extremely difficult (if possible at all) to adjust exper- 
imentally. The other possibility is to create traps which 
are almost homogenous in the middle and which have an 
appreciable trapping potential only close to the bound- 
aries. This can be studied theoretically by considering 
traps with higher powers of the trapping potentials. As 
shown below, already non-interacting systems make clear 
that a sharp Fermi edge is missing in confined systems. 

In Fig. IHIa) we show the density profile of a system 
with 1000 sites, Nf = 840, and a trapping potential of 
the form Vioxj^ with Vioa^° = 7 x lO'^'^t. It can be seen 



that the density is almost flat all over the trap with a 
density of the order of one particle per site. Only a small 
part of the system at the borders has the variation of 
the density required for the particles to be trapped. In 
Fig. IHIb) (continuous line), we show the corresponding 
normalized momentum distribution. It can be seen that 
a kind of Fermi surface develops in the system but for 
smaller values of fc, Uk is always smooth and its value 
starts decreasing at A: = 0. In order to see how nk changes 
when an incompressible region appears in the system, we 
introduced an additional alternating potential, so that in 
this case the new Hamiltonian has the form 

H = -iE(4c.+i.+H.c.) + yio^xi° n,, 

(4) 

where Va is the strength of the alternating potential. For 
the parameters presented in Fig. ISJa), we obtain that a 
small value of Va {Va = O.li) generates a band insula- 
tor in the trap, which extends over the region with n ~ 1 
(when Va = 0). However, the formation of this band insu- 
lator is barely reflected in nfc, as can be seen in Fig.|SIb) 
(dashed line). Only when the value of Va is increased and 
the system departs from the phase transition \Va = 0.5t, 
dotted line in Fig.lSfc)], does a quantitatively appreciable 
change in Uk appear. 
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FIG. 8: Exact results for Nj = 840 noninteracting trapped 
fermions in a lattice with 1000 sites and a confining potential 
Vioa^" = 7 X lO^^^f. Density profile (a) and the normalized 
momentum distribution function (b): the continuous line cor- 
responds to (a), the dashed line is the result when an alter- 
nating potential Va = Q.lt is superposed on the system, and 
the dotted line corresponds to Va = 0.5t. 
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In general, it is expected that on increasing the sys- 
tem size, the situation in the homogeneous system is ap- 
proached. However, in the present case it is not merely 
a question of boundary conditions, but the whole sys- 
tem is inhomogcneous. Therefore, a proper scaling has 
to be defined in order to relate systems of different sizes. 
In the case of particles trapped in optical lattices when 
there is a confining potential with a power a and strength 
(Vq), a characteristic length of the system (^) is given by 
C = (^a/i)~^^", so that a characteristic density (/5) can 
be defined as p = N^a/C,. We find that this characteris- 
tic density is the one meaningful in the thermodynamic 
limit. In Fig. Owe show three systems in which the to- 
tal number of particles and the curvature of the confin- 
ing potential Vio were changed, keeping p constant. We 



around a density 0.85. In the case of the momentum dis- 
tribution function [Fig.lSJb)], it is possible to see that this 
quantity also scales very well so that almost no changes 
occur in the momentum distribution function when the 
occupied system size is increased, implying that in the 
thermodynamic limit the behavior of is different from 
the one in the homogeneous system. An expanded view 
of the region with around 0.87 is introduced to better 
see the scale of the differences in this region. In Sec. V 
we show that the characteristic density defined here is 
also the meaningful quantity to define the phase diagram 
of the system. 

IV. MEAN-FIELD APPROXIMATION 
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The mean-field (MF) approach has been very useful in 
the study of dilute bosonic gases confined in harmonic 
potentials |^. This theory for the order parameter asso- 
ciated with the condensate of weakly interacting bosons 
for inhomogcneous systems takes the form of the Gross- 
Pitaevskii theory. When the bosonic condensate is loaded 
in an optical lattice, it is possible to go beyond the weakly 
interacting regime and reach the strongly correlated limit 
for which a superfluid-Mott-insulator transition occurs. 
Also in this limit a MF study was done and the re- 
sults were compared with exact diagonalization results 
for very small systems, reporting a qualitatively good 
agreement between both methods [8|. 

In this section, we compare a MF approximation with 
QMC results for the system under consideration. It will 
be shown that the MF approach not only violates the 
Mermin- Wagner theorem, leading to long-range antifer- 
romagnetic order in one dimension, as expected, but also 
introduces spurious structures in the density profiles. In 
order to obtain the MF Hamiltonian, we rewrite the Hub- 
bard Hamiltonian in Eq. in the following form (up to 
a constant shift in the chemical potential): 



FIG. 9: Exact results for noninteracting trapped fermions in 
a confining potential Vio. (a) Density profiles, (b) normalized 
momentum distribution function. Dashed line corresponds to 



Vioa^ 



6 X 10" 



,10 



428 (iV 



500), continuous line 
corresponds to Vioa'" = 7x 10"^^t, Nf = 840 (A^ ~ 1000) and 
dotted line corresponds to Vioa^° = 1 x lO'^^t, Nf = 1620 
(A'' ~ 2000). In the density profiles the positions are given in 
units of the characteristic length 



measured the positions in the trap in units of the char- 
acteristic length These systems have occupied regions 
{n > 0) with very different sizes, of the order of 500 lat- 
tice sites for the dashed line, of the order of 1000 sites for 
the continuous line, and of the order of 2000 lattice sites 
for the dotted hue. Fig.|5{a) shows that the density pro- 
files scale perfectly when the curvature of the confining 
potential and the number of fermions are changed in the 
system, so that in order to show the changes in the den- 
sity profile, we introduced an inset that expands a region 



H 



(4 
(1-A) 



H.c.) 



(5) 



where the density and the magnetization in each site are 
denoted by = fiii+nii and fci = hi^—hii, respectively. 
The fluctuations of the density and magnetization are 
given by (5„, = fii - {n^) and (5^, = /tj - (/i,;), respectively, 
and A is an arbitrary parameter that was introduced in 
order to allow for the most general variation in parameter 
space. The following relations for fermionic operators 
were used: 



1 



(6) 
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Neglecting the terms containing the square of the den- 
sity and magnetization fluctuations in Eq. (jSJ, the fol- 
lowing MF Hamiltonian 



Hmf = -t"^{clc^+la+ll.C.) + VJ2 



+ ^[(1- A) (n,)+A(M,)]n,4 
(Ij-A) A 2 



(7) 



is obtained. Due to the inhomogencity of the trapped 
system, an unrestricted Hartree-Fock scheme is used to 
determine the local densities {hi) and local magnetic 
moments {fli). Given a MF ground state \^mf), the 
minimum energy E = {"^ mf\H\'^ mf) (not the MF one 
Emf = mfWmfY^ mf)) is reached when A — 0.5, 
so that all the results that follow were obtained for this 
value of A. 

We first recall some of the discrepancies between MF 
approximations and the known valid facts for the ID ho- 
mogenous system {V = 0) as follows. 

(i) At half filling it is known that the Hubbard model 
exhibits a Mott-insulating phase with a charge gap while 
the spin sector remains gapless, so that the density- 
density correlations decay exponentially and the spin- 
spin correlations decay as a power law. Within the MF 
approximation, there is a band insulator in the system at 
half filling so that both the density-density and spin-spin 
correlations decay exponentially. The MF value of the 
charge gap is an overestimation of the real one, as can be 
seen in Fig. ^| where we present the MF results for the 
density-density correlations as function of the distance, 
together with the QMC results. The slopes of the curves 



5C 




FIG. 10: (Color online) Density-density correlations for a ho- 
mogeneous system with A'^ = 102 at half filling for U — 6t. 
MF (A) and QMC (v) solutions. 

are proportional to the charge gap and the MF slope is 
approximately twice the QMC one (the results were ob- 
tained for a system with 102 sites and U = 6t). The cor- 
relation length of X is the inverse of these slopes. Finally, 
at half filling the MF theory leads to an antiferromagnetic 
state, while in the Hubbard model the magnetization is 
always zero although quasi-long-range antiferromagnetic 
correlations appear in the system. 



(ii) At any noncommensuratc filling, the Hubbard 
model in ID describes a metal (Luttinger liquid), so that 
no gap appears in the charge and spin excitations, and 
there is a 2fci? modulation in the spin-spin correlation 
function that leads to the well-known 2fci? singularity. 
Within the MF approximation, there is always a band 
insulator at any noncommensuratc filling, with a gap 
that decreases when the system departs from half fill- 
ing. The appearance of this gap for any density is due to 
the perfect nesting present in one dimension. The insu- 
lating nature of these solutions can be seen in the global 
compressibility of the system that is always zero, in the 
behavior of the density-density and the spin-spin corre- 
lations which decay exponentially, and in the momentum 
distribution function where there is no Fermi surface, as 
shown in Fig. [TT] for iV = 102, Nj = 66, and U = 6t. 
Within this MF approach, there is a 2kp modulation of 




FIG. 11: (Color online) MF result for the momentum dis- 
tribution function of a homogeneous system with A'' — 102, 
Nf = 66, and U = 6t. 

the magnetization only in the z component [the SU (2) 
symmetry was broken] , which leads to a divergence of the 
Fourier transform of (S^S^) at fc = 2kF- In Fig. [T^ we 
compare the MF rcsuh (a) for (5'^5^)fe with the QMC 
one (b), for N = 102, Nf = 66, and U = 6t, where it 
can be seen that the 2kp peak is one order of magnitude 
bigger in the MF case compared to the QMC case. This 
is due to the existence of the magnetization in the MF 
solution (the values at the peaks will only diverge in the 
thermodynamic limit). 




Y ' ' 


0.4 r 


(a) 


0.3 - 






0.2 - 




0.1 - 




3II/4 



FIG. 12: (Color online) MF (a) and QMC (b) results for 
{S'^S^)fe in a homogeneous system with N = 102, Nf = 66, 
and U = 6t. 
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We consider next the confined inhomogeneous case. In 
Fig Jiaf a) we show the density profiles obtained for a trap 
with N = 100, Nf = 70, and Va'^ = 0.0025t (Hke the one 
presented in Fig. |3J) for different values of the on-site re- 
pulsion U/t = 2, 3, 4, 6 [for more details see Figs. lTsT eV 
I15r h)]. It can be seen that the MF solutions have also 
density profiles in which "metallic" (n ^ 1) and insulat- 
ing {n — 1) phases coexist. However, the MF insulating 
plateaus with n ~ 1 [U/t = i in Fig. llBr a)] appear for 
smaller values of U than the ones required in the Hubbard 
model for the formation of the Mott-insulating plateaus 
[U/t = 6 in Fig.|2Ia)]. In the "metallic" phases it is pos- 
sible to see that the charge density shows rapid spatial 
variations that are very large for n > l,U/t = 4 [see also 
Fig. I15f f^]. These density variations in the "metallic" 
phases are reflected in the variance profiles [Fig. Iiar b)]. 
The values of the variance in the plateaus are the same 
as the ones obtained in the homogeneous MF case for the 
same value of U when n = 1 [Fig. El (b)]. This feature 
was discussed in Sec. II for the QMC solutions. 



1.4 




FIG. 13: (Color online) MF profiles for a trap with = 
100, Va^ = 0.0025f, and Nf = 70, the on-site repulsions are 
U/t = 2 (A), 3 (v), 4 (O), and 6 (C;,). (a) Local density, (b) 
variance of the local density. The dashed lines in (b) are the 
MF values of the variance in the n = 1 homogeneous system 
for U/t = 2, 3, 4, 6 (from top to bottom). 

Since the MF solution leads to an insulating state for 
incommensurate fillings in the homogeneous case, we an- 



alyze in the following more carefully the regions with 
7T, 7^ 1 in the presence of the trap. We find that the local 
compressibility always decays exponentially as a function 
of £ over the entire system, although the exponents are 
different depending on the density of the point analyzed. 
Some additional modulation appears in when n ^ 1, 
as shown in Fig. 1141 where the local compressibility is 
displayed as a function of £ for n = 1 and n ~ 0.66 in 
the trapped and homogeneous cases. Therefore, although 
the MF approximation leads to a density profile similar 
to the one obtained in QMC as long as ti < 1, it gives 
a qualitative wrong description of the character of those 
regions. 



1 




2 4 fi 6 8 10 



FIG. 14: (Color online) MF local compressibility as a func- 
tion of £ for the homogeneous and trapped (Va^ — 0.0025f, 
Nf = 70) systems for n = 1 and n ~ 0.66, U = 6t. (v) 
n — 1 homogeneous, (Q) ni=o — 1 trapped, (A) n ~ 0.66 
homogeneous, and (O) Ui^so ^ 0.66 trapped. 

The MF results for the normalized momentum distri- 
bution function are presented in Figs. [TsTa^- tTsT d) with 
their corresponding density profiles [Figs. ll5r e)- [T5r h)]. 
for the same parameter values as in Fig. 1131 It can be 
seen that the normalized momentum profiles are similar 
to the ones obtained with QMC when the density pro- 
files are similar, so that the MF results are very similar 
to the QMC ones for this averaged quantity although the 
physical situation is very different. (Within MF there is 
always an insulator in the system and within QMC there 
are metallic and insulating phases coexisting.) In Figs. 
I15f e)- tl3r h). we also show the (S^) component of the spin 
on each site of the trap [S^ = {hi^ — fiii) /2]. When the 
density is around n = 1, it can be seen that antiferro- 
magnetic order appears, and that a different modulation 
exists for n ^ 1. We should mention that we have also 
found some MF solutions for the trapped system where 
there was formation of spin domain walls in the insulat- 
ing regions with n = I and there the variance was not 
constant (a plateau in the density does not correspond to 
a plateau in the variance), so that care should be taken 
with the mean field solutions. 
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FIG. 15: (Color online) MF normalized momentum distribution function (Q) and their corresponding density (A) and 

(5"> (V) profiles (e)-(h) for U/t = 2 (a),(e), 3 (b),(f), 4 (c),(g), 6 (d),(h) and = 100, Nf = 70, Va^ = 0.0025*. 



V. PHASE DIAGRAM 

In the present section, we study the phase diagram for 
fermions confined in harmonic traps with an underlying 
lattice. 

As it can be inferred from Fig. ^ the phases present in 
the system are very sensitive to the values of the param- 
eters in the model, i.e., to the curvature of the parabolic 
potential (V), to the number of fermions present in the 
trap (Nf) and to the strength of the on-site repulsion 



(f/). In order to be able to relate systems with different 
values of these parameters and different sizes, we use the 
characteristic density {p) defined in Sec. Ill for the defini- 
tion of the phase diagram. This quantity has been shown 
to be a meaningful quantity in the thermodynamic limit 
since density profiles (as function of x/C) and normalized 
momentum distributions do not change when jj is kept 
constant and the filling in the system (or the occupied 
system size) is increased. For a harmonic potential, the 

characteristic length is given by C = {V/ty^ and the 
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1/2 

characteristic density is then p = Nfa{V/t) ' . In Fig. 
1161 we show the phase diagrams of two systems with dif- 
ferent strengths of the confining potential {Va^ = 0.006t 
and Va'^ = 0.002i) and different sizes {N = 100 and 
N — 150, respectively). This shows that the character- 
istic density is a meaningful quantity to characterize a 
generic phase diagram. It allows to compare different 
systems, and, hence, to relate the results of numerical 
simulations with larger experimental sizes (although cur- 
rently the linear dimension of the experiments is ^ 65 a). 
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FIG. 16: (Color online) Phase diagram for systems with 
Va'^ = 0.0064, N = 100 (v) and Va^ = 0.002t, iV = 150 
(O)- The different phases are explained in the text. 

The phases present in the phase diagram in Fig. ll6l can 
be described as follows: (A) a pure metal without insu- 
lating regions, (B) a Mott insulator in the middle of the 
trap always surrounded by a metal, (C) a metallic region 
with n > 1 in the center of the Mott-insulating phase, 
that is surrounded by a metal with n < 1, (D) an insula- 
tor with n = 2 in the middle of the trap surrounded only 
by a metallic phase, and finally (E) a "band insulator" in 
the middle of the trap surrounded by metal with n > 1, 
and these two phases inside a Mott insulator that as al- 
ways is surrounded by a metal with n < 1. As it can be 
seen, this phase diagram is more complex than the one 
for the ID homogeneous Hubbard model. 

There are some features that we find interesting in this 
phase diagram: (i) For all the values of U that we have 
simulated, we see that the lower boundary of phase B 
with the phase A has a constant value of the characteris- 
tic density, (ii) We have also seen that the upper bound- 
ary of phase B with phases A and C and the boundary 
of phases D and E with A and C, respectively, are linear 
within our errors, (iii) Finally, there is another character- 
istic of phase B, and of the metallic phase A that is below 
phase B, that we find intriguing and could be related to 
point (i). As can be seen in Fig. |3{a) for U/t ~ 6 and 8, 
once the Mott-insulating phase is formed in the middle of 
the trap, a further increase of the on-site repulsion leaves 
the density profile almost unchanged. (The two men- 
tioned profiles arc one on top of the other in that figure.) 



However, if we look at the variance of these profiles it 
decreases with increasing U, which means that the dou- 
ble occupancy = [A — n (1 — n)] /2 is decreasing 
in the system when U is increased, without a redistribu- 
tion of the density. Less expected for us is that this also 
occurs for the metallic phase A that is below phase B. 
In Fig. llTr a) we show four density profiles and in Fig. 
I17f b'l their variances, for a trap with Va^ — 0.006t and 
Nf = 30 when the on-site repulsion has values U/t = 2, i, 
6, and 8. If we look at the phase diagram, the point with 
U/t = 2 is not below phase B whereas the other three 
values are. In Fig. llTf a) we can see that as the value of 
U/t is increased from 2 to 4 there is a visible change in 
the density profile, but after a further increase the profile 
remains almost the same for the other values of U. The 
inset shows a magnification of the top of the profiles for 
a more detailed view. Fig. llTf b') shows that although the 
density stays almost constant for U/t = 4, 6, and 8, the 
variance decreases. In the other phases (C, D, and E), an 
increase of the on-site repulsion always changes the local 
densities pushing particles to the edges of the trap. 




FIG. 17: (Color online) Profiles for a trap with Va^ = 0.006t 
and Nf = 30, the on-site repulsions are U/t = 2 (A), 4 (v)i 
6 (O), and 8 (O)- (a) Local density, (b) variance of the local 
density. 

Finally, in order to sec, whether the features of the 
phase diagram discussed previously arc particular of a 
harmonic confining potential, we have also performed an 
analogous study for a quartic confining potential. In Fig. 
I18l we show how the density profiles evolve in a trap with 
a quartic confining potential [V^^a^ = 3 x IQ-^t) when 
the total filling is increased from Nf ~ 20 to 140. It 
can be seen that the shape of the metallic regions is fiat- 
ter than in the harmonic case. However, all the features 
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discussed previously for the parabolic confining poten- 
tial are also present in the quartic case. Local phases 
appear in the system, the Mott-insulating plateaus with 
n = 1 have values of the variance equal to the ones in 
the homogeneous case for the same value of U, the local 
compressibility vanishes in the Mott-insulating regions 
and the spin and momentum distribution exhibit a be- 
havior similar to that of the parabolic case. In addition, 




FIG. 18: Evolution of the local density in a quartic confining 
potential as a function of the position in the trap and in- 
creasing total number of fermions. The parameters involved 
are iV = 100, U = 5t, and l^a" = 3 x IQ-^t. 

the form of the phase diagram for the quartic confin- 
ing potential is similar to that of the harmonic trap. In 
Fig. ^1 we show the phase diagram for a system with 
N = 100 and Viu'^ = 3 x lO^^t, where the characteris- 
tic density is given by p = Nja {V4/tY^^ . We have also 
calculated some points of the phase diagram for another 
system with V4a^ = lO^^t and N = 100 in order to check 
that the scaling relation also works properly in this case. 
The phases arc labeled in the same way as in Fig. 1161 up 
to the largest value of U that we have studied we do not 
obtain the phase E which for the quartic confining case 
seems to be moved further towards the strong coupling 
regime. Figure [T^ shows that the phase diagram in this 
case is very similar to that of the parabolic case, so that 
all the features discussed previously apply to the quartic 
confining potential. The form of the phase diagram in 
Fig. 1161 seems to be of general applicability for any con- 
fining potential when the proper characteristic density is 
considered. 



VI. CONCLUSIONS 

We have studied the ground-state properties of the ID 
fermionic Hubbard model in harmonic traps with an un- 
derlying lattice, using QMC simulations and a MF ap- 
proach. 

Due to the inhomogencous density distribution, in gen- 
eral metallic and insulating regions coexist in the trap. 
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FIG. 19: (Color online) Phase diagram for systems with 
Via"^ = 3 X lO'^t, N = 100 (v) and Via"* = lO'^'t, N = 100 
(O)- The phases are labeled in the same way as in Fig. 1161 



Therefore, local quantities have to be used to characterize 
the system. We considered first the variance of the local 
density, which, when a plateau develops in the density 
profile with n = 1, equals the value in the homogeneous 
system in the Mott-insulating phase. However, it is not 
an unambiguous quantity, since it can be larger in the 
Mott-insulating phases than in metallic regions with den- 
sities close to n or n ~ 2. Therefore, we have defined 
a local-order parameter (local compressibility) that van- 
ishes in the insulating regions and is finite in the metallic 
ones. The local compressibility gives a measure of the 
local change in density due to a constant shift of the po- 
tential over finite distances larger than the correlation 
length of the density-density correlation function for the 
Mott-insulating phase in the unconfincd system. As ex- 
pected, antiferromagnetic correlations are present in the 
Mott-insulating phases and they remain in the shoulders 
of the density profile for n ^ 1. 

As opposed to homogeneous or periodic systems, where 
the appearance of a gap can in general be seen by the dis- 
appearance of the Fermi edge in the momentum distribu- 
tion function, it is easily seen that due to the presence of 
a confining potential, such a feature is much less evident. 
Although increasing the power of the confining poten- 
tial sharpens the features connected with the Fermi edge, 
it remains qualitatively different from the homogeneous 
case. Increasing the system size with a proper definition 
of the thermodynamic limit, does not change at all the 
behavior of the density and momentum distribution func- 
tions in the trapped case, such that, in fermionic systems 
it does not seem to be the appropriate quantity to look 
for evidence of gaps opening in the system. 

The MF study has shown that within this approach 
the trapped system is an insulator for all the densities 
in the trap, a qualitatively wrong picture in contradic- 
tion with the QMC results. In addition, the metallic 
regions showed large oscillations in the density that are 
not present in the real solutions. However, it is possible 
to obtain some qualitative information from MF, such as 
the coexistence of n 7^ 1 regions with n = 1 plateaus, the 
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shape of the momentum distribution function and the ex- 
istence of antiferromagnetic correlations in the insulating 
regions with n — I. 

Finally, we have determined a generic form for the 
phase diagram that allows to compare systems with dif- 
ferent values of all the parameters involved in the model. 
It can be used to predict the phases that will be present 
in future experimental results. The phase diagram also 
reveals interesting features such as reentrant behavior 
in some phases when some parameters are changed and 
phase boundaries with linear forms. A similar form of 
the phase diagram has been also found for a quartic con- 
fining potential, such that the results obtained here are 
generally applicable to cases beyond a perfect harmonic 
potential. 
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